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Abstract 

A robust, fast and accurate method for solving the Colebrook- hke 
equations is presented. The algorithm is efficient for the whole range of 
parameters involved in the Colebrook equation. The computations are not 
more demanding than simplified approximations, but they are much more 
accurate. The algorithm is also faster and more robust than the Colebrook 
solution expressed in term of the Lambert W-function. Matlab® and 
FORTRAN codes are provided. 

1 Introduction 

Turbulent fluid flows in pipes and open channels play an important role in 
hydraulics, chemical engineering, transportation of hydrocarbons, for example. 
These flows induce a significant loss of energy depending on the flow regime 
and the friction on the rigid boundaries. It is thus important to estimate the 
dissipation due to turbulence and wall friction. 

The dissipation models involve a friction coefficient depending on the flow 
regime (via a Reynolds number) and on the geometry of the pipe or the channel 
(via an equivalent sand roughness parameter) . This friction factor if often given 
by the well-known Colebrook-White equation, or very similar equations. 

The Colebrook-White equation estimates the (dimcnsionless) Darcy-Weis- 
bach friction factor A for fluid flows in flUed pipes. In its most classical form, 
the Colebrook-White equation is 

1 „ , f K 2.51 1 \ 

where R — UD/v is a (dimcnsionless) Reynolds number and K ~ e/D \s a. 
relative (dimcnsionless) pipe roughness {U the fluid mean velocity in the pipe, 
D the pipe hydraulic diameter, v the fluid viscosity and e the pipe absolute 
roughness height). There exist several variants of the Colebrook equation, e.g. 
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These variants ean be recast into the form (1) with small changes in the numer- 
ical constants 2.51 and 3.7. Indeed, the latter numbers being obtained fitting 
experimental data, they are known with limited accuracy. Thus, the formulae 
(2) and (3) are not fundamentally different from (1). Similarly, there are vari- 
ants of the Colebrook equations for open channels, which are very similar to (1). 
Thus, we shall focus on the formula (1), but it is trivial to adapt the resolution 
procedure introduced here to all variants, as demonstrated in this paper. 

The Colebrook equation is transcendent and thus cannot be solved in terms 
of elementary functions. Some explicit approximate solutions have then been 
proposed [6, 10, 12]. For instance, the well-known Haaland formula [6] reads 



Haaland's approximation is explicit but is not as simple as it may look. Indeed, 

this approximation involves one logarithm only, but also a non-integer power. 
The computation of the latter requires the evaluation of one exponential and 
one logarithm, since it is generally computed via the relation 



where 'In' is the natural (Napierian) logarithm. Hence, the overall evaluation 
of (4) requires the computation of three transcendant functions (exponentials 
and logarithms). We present in this paper much more accurate approxima- 
tions requiring the evaluation of only two or three logarithms, plus some trivial 
operations (+, — , x, -^). 

Only quite recently, it was noticed that the Colebrook- White equation (1) 
can be solved in closed form [8] using the long existing Lambert W- function [3]. 
However, when the Reynolds number is large, this exact solution in term of the 
Lambert function is not convenient for numerical computations due to overflow 
errors [11]. To overcome this problem, Sonnad and Goudar [11, 12] proposed 
to combine several approximations depending on the Reynolds number. These 
approaches are somewhat involve and it is actually possible to develop a simpler 
and more efficient strategy, as we demonstrate in this paper. 

A fast, accurate and robust resolution of the Colebrook equation is, in par- 
ticular, necessary for scientific intensive computations. For instance, numerical 
simulations of pipe flows require the computation of the friction coefficient at 
each grid point and for each time step. For long term simulations of long pipes, 
the Colebrook equation must therefore be solved a huge number of times and 
hence a fast algorithm is required. An example of such demanding code is the 
program OLGA [1] which is widely used in the oil industry. 

Although the Colebrook formula itself is not very accurate, its accurate 
resolution is nonetheless an issue for numerical simulations because a too crude 
resolution may affect the repeatability of the simulations. Robustness is also 
important since one understandably wants an algorithm capable of dealing with 
all the possible values of the physical parameters involved in the model. The 
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exp(l.ll X ln(.x)). 
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method described in the present paper was developed to address aU this issues. 
It is also very simple so it can be used for simple applications as well. The 
method proposed here aims at giving a definitive answer to the problem of 
solving numerically the Colebrook-like equations. 

The paper is organized as follow. In section 2, a general Colebrook-like 
equation and its solution in term of the Lambert W-function are presented. For 
the sake of completeness, the Lambert function is briefly described in section 
3, as well as a standard algorithm used for its computation. A severe draw- 
back of using the Lambert function for solving the Colebrook equation is also 
pointed out. To overcome this problem, a new function is introduced in section 
4 and an improved new numerical procedure is described. Though this func- 
tion introduces a big improvement for the computation of the friction factor, 
it is still not fully satisfactory for solving the Colebrook equation. The reasons 
are explained in the section 5, where a modified function is derived to address 
the issue. The modified function is subsequently used in section 6 to solve the 
Colebrook equation efficiently. The accuracy and speed of the new algorithm 
is tested and compared with Haaland's approximation. For testing the method 
and for intensive practical applications, Matlab© and FORTRAN implemen- 
tations of the algorithm are provided in the appendices. The algorithm is so 
simple that it can easily be implemented in any other language and modified to 
be adapted to any variants of the Colebrook equation. 



2 Generic Colebrook equation and its solution 

We consider here a generic Colebrook-like equation as 

= CO - ciln(c. + ^), (5) 

where the Ci are given constants such that C1C3 > 0. The classical Colebrook- 
White formula (1) is obviously obtained as a special case of (5) with co = 0, 
ci = 2/ In 10, C2 = and C3 = 2.51/i?. 

The equation (5) has the exact analytical solution 



7t ^ 



Wicxpl £^l + ^_ln(ciC3)^ ^ 



(6) 



Cl C1C3 J J C1C3 

which is real if C1C3 > and where W is the principal branch of the Lambert 
function, often denoted Wq [3]. In this paper, only the principal branch of the 
Lambert function is considered because the other branches correspond to non- 
physical solutions of the Colebrook equations, so the simplified notation W is 
not ambiguous. 



3 Brief introduction to the Lambert W-function 

For the sake of completeness, we briefly introduce the Lambert function and its 
practical computation. Much more details can be found in [3, 4] . 
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The Lambert W-function solves the equation 



y exp(j/) = X y = W(a;), (7) 

where, here, x is real — more precisely a; > — exp(— 1) — and W(0) = 0. The 

Lambert function cannot be expressed in terms of elementary functions. An 
efficient algorithm for its computation is based on Halley's iterations [3] 

y.^^ ^ y. VJ eXp(?/j-) - X ^ 

^"^^ ^ (% + 1) cxp(yj) - \{yj + 2) {yj exp(yj) - x)/{yj + 1) ' 

provided an initial guess y^. Halley's method is cubic (c.f. Appendix A), mean- 
ing that the number of exact digits is (roughly) multiplied by three after each 
iteration. Today, programs for computing the Lambert function arc easily found. 
For instance, an efficient implementation in Matlab® (including complex ar- 
gument and all the branches) is freely available [5] . 

The Taylor expansion around a; = of the Lambert function is 

oo / \n— 1 

W(^) = E I kl<exp(-l). (9) 

n! 

n=l 

This expansion is of little interest to solve the Colebrook equation because, in 
this context, the corresponding variable x is necessarily large (a; 1). It is thus 
more relevant to consider the asymptotic expansion 

W(x) ^ ln(x) — ln(ln(x)) as x — > oo. (10) 

This expansion reveals that W behaves logarithmically for large x, while we must 
compute W(exp(a;)) to solve the Colebrook equation, c.f. relation (6). For our 
applications, x is large and exp(a;) is therefore necessarily huge, to an extend that 
the computation of cxp(a;) cannot be achieved due to overflow. Even when the 
intermediate computations can be done, the result can be very inaccurate due 
to large round-off errors. Therefore, the resolution of the Colebrook equation 
via the Lambert function [8] is not efficient for the whole range of parameter of 
practical interest [11]. 

4 The a;-function 

To overcome the numerical difficulties related to the Lambert W-function, when 
used for solving the Colebrook-White equation, we introduce here a new func- 
tion: the w-function. 

The w-function is defined such that it solves the equation 

y + ln(y) = x =^ y = uj{x), (11) 

where we consider only real x. The cj-function is related to the W-function as 

oj{x) = W(exp(x)). 
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Note that the Lambert W-function is also sometimes caUed the Omega function, 
that should not be confused with the w-function defined here, where we follow 
the notation used in [9]. In terms of the w-function, the solution of (5) is of 
course 

1 



7^ = 



(12) 
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For large arguments uj{x) behaves like x, i.e. we have the asymptotic behavior 

ijli{x) ~ X — ln(x) as x ^ oo, (13) 

which is an interesting feature for the application considered in this paper. 

As noted by Cor less et al. [4], the equation (11) is in some ways nicer than 
(7). In particular, its derivatives (with respect of y) are simpler, leading thus to 
algebraically simpler formulae for its numerical resolution. An efficient iterative 
quartic scheme (c.f. Appendix A) is thus 



with 



Vj + MVj) - X 1 

" TT^^ ' ^0 = - - 5- 



The computationally costless initial guess [yo = x—^) was obtained considering 
the asymptotic behavior (10), minus an empirically found correction (the term 
— ^) to improve somewhat the accuracy of yo for small x without affecting the 
accuracy for large x. The relative error Cj of the j-th iteration, i.e. 



i{x) = 



yjjx) - uj{x) 

LO{x) 



is displayed on the figure 1 for 1 ^ a; ^ 10^ and j = 0, 1, 2. (The accuracy of (14) 
were measured using the arbitrary precision capability of Mathematica® .) We 
can see that with j = 2we have already reached the maximum accuracy possible 
when computing in double precision, since max(e2) « 4 x 10""'^^ for x G [l;oo[. 
We note that the relative error continues to decay monotonically as x increases 
(even for x > 10^) and that there are no overflow problems when computing 
yj even for very large x (i.e. x ^ 10^). We note also that for x ^ 5700 the 
machine double precision is obtained after one iteration only. 

The scheme (14) is quartic, meaning that the number of exact digits is mul- 
tiplied by four after each iteration (c.f. Appendix A). Hence, starting with 
an initial guess with one correct digit, four digits are exact after one iteration 
and sixteen after two iterations. That is to say that the machine precision (if 
working in double precision) is achieved after two iterations only (Fig. 1). More- 
over, the scheme (14) has a comparable algebraic complexity per iteration than 
the scheme (8), i.e. the computational times per iteration are almost identi- 
cal. However, the iterative quartic scheme (14) converges faster than the cubic 
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one (8), and there are no overflow problems as they appear when computing 
W(exp(a;)) for large x. This algorithm could therefore be used to compute the 
solution of the Colebrook-White equation (1), but we will use instead an even 
better one defined in the next section. We note in passing that the iterations 
(14) are also efficient for computing the w-function for any complex x, provided 
some changes in the initial guess yo depending on x. 

Remarks: 

ir With a more accurate initial guess yo) such asyo = x — ln(a;), the desired 
accuracy may be obtained with fewer iterations. However, the computation of 
such an improved initial guess requires the evaluation of transcendent functions. 
Thus, it cannot be significantly faster than the evaluation of y\ with (14) from 
the simplest giiess = x — ^, and most likely less accurate. 

ii- Higher-order iterations are generally more involved per iteration than 
the low-order ones. Higher-order iterations are thus interesting if the number 
of iterations is sufficiently reduced so that the total computation is faster to 
achieve the desired accuracy. This is precisely the case here. 

iiir Intensive tests have convinced us that the choice of the simplest initial 
guess 7/0 = ^ ^ together with the quartic iterations (14) is probably the best 
possible scheme for computing the w-function in the interval x € [1; oo[, at least 
when working in double precision. If improvements can be found, they are thus 
most likely very minor in terms of both robustness, speed and accuracy. 

5 The tu-function 

Solving the Colcbrook equation via the cj-function is a big improvement com- 
pared to its solution in term of the Lambert W-function. One can check that the 
numerical resolution of the Colebrook equation via the algorithm (14) is indeed 
very efficient when K = Q, even for very large R. However, when K > Q the 
scheme (14) is not so effective for large R, meaning that not all the numerical 
shortcomings have been addressed introducing the w-function. The cause for 
these numerical problems can be explained as follow. 

The solution of the Colebrook equation requires the computation of an ex- 
pression like u){x\ + x^) — x\ where x\ ^ X2 when R is large and K ^ (but 
xi = ii K = 0), sec the relation (19) below. Assuming X2 oc ln(a;i), as is the 
case here, the asymptotic expansion as Xi — > oo, i.e. 

w{xi+X2) — xi ~ (xi + a;2 — ln(a;i)) — xi = X2 — In(xi), 

exhibits the source of the numerical problems. Indeed, when K > and R is 
large, we have Xi ^ X2 and Xi ^ ln(a;i). Therefore \x2 — ln(a;i)|/a;i can be 
smaller than the accuracy used in the computation and we thus obtain numeri- 
cally X1+X2 — In(xi) ss xi due to round-off errors. Hence uj{xi +X2) — .xi ss is 
computed instead of lo{xi + X2) — Xi « X2 — ln(a;i). To overcome this problem 
we introduce yet another function: the ro-function. 
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Introducing the change of variable y = z + xi into the equation (11), the 
tu-function is defined such that it solves the equation 



z + ln(xi+^) = X2 => z = 'a7{xi\x2), (15) 
where the real. The raj-function is related to the lj- and W-functions as 

Tu{xi\x2) = ti;(a;i+a;2) - xi = W(exp(xi + 0:2)) - xi. 
In terms of the ro-function, the solution of (5) is obviously 



1 f C2 



^ - ln(ciC3) ) . (16) 

Cl 



The ro-function is nothing more than the w-function shifted by the quantity Xi . 
This is a very minor analytic modification but this is a numerical significant 
improvement when xi is large. 

An efficient numerical algorithm for computing the oj-function is directly 
derived from the scheme (14) used for the w-function. We thus obtain at once 

- . (l + a;i+Zj- + ^ej-)ej-(xi+Zj-) 

- - {l + x,+z, + e, + ^ef) ' ^ ^''^ 



with 



_ Zj + ln{xi + Zj) - X2 1 

= 1 + XI + z, ' ^0 = X2- -. 



^3 

If Xi = the scheme (14) is recovered. The rate of convergence of (17) is of 
course identical to the scheme (14). Thus, the efficiency of (17) does not need 
to be re-discussed here (see section 4). 



6 Resolution of the Colebrook— White equation 



We test the new procedure with the peculiar Colebrook- White equation (1). Its 
general solution is 



1 

7x 
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18.574 
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5.02 

£R 



5.02/77 18.574 
IKR 



18.574 



5.02 



(18) 
(19) 
(20) 



where £ = In(lO) « 2.302585093. AH these analytic solutions are mathematically 
equivalent, but the relation (20) is more efficient for numerical computations if 
we use the scheme described in the previous section. 
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6.1 Numerical procedure 



The solution of the Colebrook-White equation is obtained computing the vj- 
function with 



and using the iterative scheme (17) with j = 0, 1,2. An approximation of the 
friction factor is eventually 



This way, the whole computation of \j requires the evaluation of j + l logarithms 
only,^ i.e. one logarithm per iteration. 

A Matlab® implementation of this algorithm is given in the appendix B. 
This (vectorized) code was written with clarity in mind, so that one can test 
and modify easily the program. This program is also fast, accurate and robust, 
so it can be used in real intensive applications developed in Matlab. 

A FORTRAN implementation of this algorithm is given in the appendix 
C. This program was written with speed in mind, so there are no checks of the 
input parameters. The code is clear enough that it should be easy to modify 
and to translate into any programming language. 

6.2 Accuracy 

For the range of Reynolds numbers 10^ ^ ^ 10^^ and for four relative rough- 
ness K = {0, 10~^, 10~^, 10~^}, the accuracy of A~^^^ — obtained from the 

— 1/2 

iterations (17) with j = {0, 1, 2} — and of Haaland's approximation Ajj — 
given by (4) — are compared with the exact friction coefficient A~^/^. The 
relative errors are displayed on the figure 2. 

It appears clearly that Aj is accurate to machine double-precision (at 
least) for all Reynolds numbers and for all roughnesses (in the whole range of 
physical interest, and beyond). 

— 1/2 

It also appears that Aj^ is more accurate than Haaland's approximation, 
specially for large R and K. Moreover, the computation of Ai requires the eval- 
uation of only two logarithms, so it is faster than Haaland's formula. Note that 
other explicit approximations having more or less the same accuracy as Haa- 

— 1 /2 . ■ ■ 

land's formula, is significantly more accurate than these approximations. 

— 1/2 

Finally, we note that Ag is a too poor approximation to be of any practical 



Testing the actual speed of an algorithm is a delicate task because the run- 
ning time depends of many factors independent of the algorithm itself (imple- 

^The numerical constant ln{10) is not counted because it can be explicitly given in the 
program and does not need to be computed each time. 



18.574' 




A, « {ll2z,f. 



interest. 



6.3 Speed 
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mentation, system, compiler, hardware, etc.), specially on multi-tasking and 
multi-users computers. In order to estimate the speed of our scheme as fairly 
as possible, the following methodology was used. 

The speeds of the computation of Ai and A2 are compared with the Haa- 
land approximation Ah- The Matlab environnement and its built-in cputime 
function is used, for simplicity. 

Two vectors of N components, with 1 ^ ^ 10"'', arc created for R and K. 
The values are chosen randomly in the intervals 10'^ ^ -R ^ 10^ and ^ K < 1. 
The computational times are measured several times, the different procedures 
being called in different orders. For each value of N, the respective timings are 
averaged and divided by the averaged time used by the Haaland approximation 
(the latter having thus a relative computational time equal to one for all N). 
The result of this test are displayed on the figure 3. (The whole procedure was 
repeated several times and the corresponding graphics were similar.) 

For small N, say A'' < 2000, the computations are so fast that the function 
cputime cannot measure the times. For larger values of A^, we can sec on the 
figure 3 that the computations of Ai are a bit faster than the Haaland formula, 
while the computations of A2 are a bit slower, in average. This is in agree- 
ment with the number of evaluations of transcendent functions needed for each 
approximations, as mentioned above. 

These relative times may vary depending on the system, hardware and soft- 
ware, but we believe that the results would not be fundamentally different from 
the ones obtained here. The important result is that the procedure presented 
in this paper is comparable, in term of speed, to simplified formulae such as 
the Haaland approximation. The new procedure being much more accurate, it 
should thus be preferred. 

7 Conclusion 

We have introduced a simple, fast, accurate and robust algorithm for solving 
the Colebrook equation. The formula used is the same for the whole range 
of the parameters. The accuracy is around machine double precision (around 
sixteen digits). The present algorithm is more efficient than the solution of 
the Colebrook equation expressed in term of the Lambert W- function and than 
simple approximations, such as the Haaland formula. 

We have also provided routines in Matlab and FORTRAN for its practical 
use. The algorithm is so simple that it can easily be implemented in any other 
language and can be adapted to any variant of the Colebrook equation. 

To derive the algorithm, we introduced two special functions: the ui- and 
tu-functions. These functions could also be useful in other contexts than the 
Colebrook equation. The efficient algorithms introduced in this paper for their 
numerical computation could then be used, perhaps with some modifications of 
the initial guesses, specially if high accuracy is needed. 
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A High-order schemes for solving a single non- 
linear equation 

Let be a single nonlinear equation f(y) = 0, where / is a sufBciently regular 
given function and y is unknown. This equation can be solved iteratively via 
the numerical scheme [7] 

(21) 

V=Vi 

where p is a non- negative integer and F^p^ denotes the p-th derivative of F with 
i?(0) = F 

The scheme (21) is of order p + 2, meaning that the number of exact digits is 
roughly multiplied by p + 2 after each iteration (when the procedure converges, 
of course). For p = and p = I, one obtains Newton's and Halley's schemes, 
respectively. The scheme (14) for solving the Colebrook equation is obtained 
with p = 2 together with the function / given by the equation (11), plus some 
elementary algebra. Intensive tests have convinced us that it is most probably 
the best choice for the problem at hand here. 



%+i = Vj + (P+1) 



(1//) 
(l//)(P+i)_ 



B MATLAB code 



The Matlab® function below is a vectorized implementation of the algorithm 
described in this paper. This code can also be freely downloaded [2]. We hope 
that the program is sufficiently documented so that one can easily test and 
modify it. 



function F = colebrook (R,K) 

% F = COLEBR]0OK(R,K) fast , accurate and robust computation of the 
% Darcy— Weisbach friction factor according to the Colebrook formula; 
% - - 

% 1 I K 2.51 I 

% = -2 * LoglO I + I 

% sqrt(F) I 3.7 R * sqrt(F) | 

% - - 

% INPUT: 

% R : Reynolds' number (should be > 2300). 

% K : Equivalent sand roughness height divided by the hydraulic 
% diameter (default K=0). 

% 

% OUTPUT: 

% F : Friction factor. 
% 

% FORMAT: 

% R, K and F are either scalars or compatible arrays. 

% 

% AOCURACY: 

% Around machine precision for all R > 3 and for all <— K, 

% i.e. in an interval exceeding all values of physical interest . 

% 

% EXAMPLE: F= colebrook ([ 3 e3 , 7 e5 , 1 elOO ], . 1 ) 

% Check for errors . 
if any(R(:)<=0) = 1, 
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error ('The Reynolds number must be 
end , 

i f nargin = 1 , 

K = 0; 
end , 

if any (K(:) <0) = 1 , 

error ('The relative sand roughness 
end , 

% Initialization. 

XI = K .* R * 0.123968186335417556; 
X2 = log(R) - 0.779397488455682028; 

% Initial guess. 
F = X2 - 0.2; 



positive (R>2000).'); 

must be non— negative .') ; 

% XI <- K * R * log (10) / 18.574. 
% X2 <- log( R * log(lO) / 5.02 ); 

% F <- X2 - 1/5; 



% First iteration. 

E = ( log(Xl+F) + F-X2 ) ./ ( 1+Xl + F ); 

F = F - (1+X1+F + 0.5*E) .* E .* (Xl+F) ./ (1+Xl+F+E. * (1+E/ 3 ) ) ; 

% Second iteration (remove the next two lines for moderate accuracy). 
E = ( log (Xl+F) + F-X2) ./ (1+Xl + F); 

F = F - (1+X1+F + 0.5*E) .* E .» (Xl+F) ./ (1+Xl+F+E. * (1 +E/ 3 ) ) ; 
% Finalized solution. 

F = 1.151292546497022842 ./ F; % F <- 0.5 * log(lO) / F; 

F = F .* F; % F <- Friction factor. 



C FORTRAN code 



The FORTRAN function below was written with maximum speed in mind, so 
some trivial arithmetic simplifications were used and there are no check for 
errors in the input parameters. 



DOUBLE PRECISION FUNCTION COLEBROOK(R,K) 

C F = COLEBROOK(R,K) computes the Darcy-Weisbach friction 

C factor according to the Colebrook— White formula. 

C 

C R : Reynold ' s number. 

C K : Roughness height divided by the hydraulic diameter . 

C F : Friction factor. 

IMPLICIT NONE 

DOUBLE PRECISION R, K, F, E, XI, X2 , T 
PARAMETER ( T = 0.333333333333333333D0 ) 

C Initialization. 

XI = K * R * 0. 123968186335417556D0 
X2 = LOG(R) - 0.779397488455682028D0 

C Initial guess . 

F = X2 - 0.2D0 

C First iteration. 

E= (LOG(Xl+F) -0.2D0) / ( 1 . DOfXl+F ) 

F = F - (1.0D(>fXl+F+0.5D0*E)*E*(Xl+F) / ( 1 . ODOfXl+F+E* ( 1 . ODafE*T) ) 

C Second iteration (if needed). 

IF ((X1+X2) .LT. (5 . 7D3)) THEN 

E = (LOG(Xl+F)+F-X2) / (l.ODOfXl+F) 

F = F - (1.0DafXl+F + 0.5D0*E)*E*(Xl+F) / ( 1 . ODOfXl+F+E* ( 1 . OD(H-E*T) ) 
ENDIF 
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C Finalized solution . 

F = 1.151292546497022842D0 / F 
(X>LEBROOK = F * F 



RETURN 
END 

Note that, depending on the FORTRAN version and on the compiler, the 
command LOG may have to be replaced by DLOG to ensure that the logarithm 
is computed with a double-precision accuracy. 
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Figure 1: Relative errors Cj of the to-function computed via the iterations (14)- 
Dotted red line: eo; Dashed blue line: ei; Solid green line: 62- 
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Figure 2: Relative errors of , computed via the iterations (17) and the 
Haaland formula (4), as functions of the Reynolds number R. Upper-left: K = 
0; Upper-right: K — ; Lower-left: K = 10~^; Lower-right: K ~ 10^^. 

Dotted red line: Aq ^ ; Dashed blue line: Xi ^ ; Solid green line: A2 ^ ; 

_ 1 

Dashed-dotted black line Ajj^ (Haaland's approximation). 
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Figure 3: Computational times of Xi (dashed blue line) and A2 (solid green line) 
with respect of the Haaland approximation Xh (dashed- dotted black line). 
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